Comparison: bayesQR vs tauBayesW Methods

tauBayesW Package

2025-09-10

Introduction

This vignette compares the performance of bayesQR with the three MCMC methods implemented in tauBayesW using the classic Prostate cancer dataset. We’ll estimate the 75th quantile (τ = 0.75) and compare the coefficient estimates across all methods.

The tauBayesW package implements three MCMC approaches:

tauBayesW Usage Examples

Before diving into the comparison, let’s see how to use tauBayesW methods:

Creating Priors

# Create a prior for p=5 coefficients
my_prior <- tauBayesW::prior(
  p = 5,                                    # Number of parameters
  type = "univariate",                      # Prior type
  beta_mean = rep(0, 5),                   # Prior mean for coefficients
  beta_cov = diag(1000, 5),               # Prior covariance matrix
  sigma_shape = 0.001,                     # Shape parameter for sigma (ALD only)
  sigma_rate = 0.001,                      # Rate parameter for sigma (ALD only)
  names = c("Int", "x1", "x2", "x3", "x4") # Parameter names
)
print(my_prior)

Fitting Models with Different Methods

# Example data
data(mtcars)

# Method 1: ALD (Asymmetric Laplace Distribution)
fit_ald <- bqr.svy(
  mpg ~ wt + hp + cyl,
  data = mtcars,
  quantile = 0.5,
  method = "ald",
  niter = 10000,
  burnin = 5000,
  thin = 5,
  print_progress = 2500  # Print every 2500 iterations
)

# Method 2: Score Likelihood
fit_score <- bqr.svy(
  mpg ~ wt + hp + cyl,
  data = mtcars,
  quantile = 0.5,
  method = "score",
  niter = 10000,
  burnin = 5000,
  thin = 5,
  print_progress = 2500
)

# Method 3: Approximate Method
fit_approx <- bqr.svy(
  mpg ~ wt + hp + cyl,
  data = mtcars,
  quantile = 0.5,
  method = "approximate",
  niter = 10000,
  burnin = 5000,
  thin = 5,
  print_progress = 2500
)

Model Output: Summary and Print Methods

# Print method - shows basic model information
print(fit_ald)

# Summary method - detailed convergence diagnostics
summary(fit_ald)

# Extract coefficient estimates
coef_estimates <- fit_ald$beta
print(coef_estimates)

# Access posterior draws
posterior_draws <- fit_ald$draws
head(posterior_draws)

# Check acceptance rate (if available)
if (!is.null(fit_ald$accept_rate)) {
  cat("Acceptance rate:", fit_ald$accept_rate, "\n")
}

Data Preparation

library(tauBayesW)
library(bayesQR)

# Load Prostate dataset
data(Prostate, package = "bayesQR")

# Center covariates (subtract mean, don't scale)
covars <- Prostate[, colnames(Prostate) != "lpsa"]
covars_centered <- as.data.frame(scale(covars, center = TRUE, scale = FALSE))

# Combined dataset with centered covariates
Prostate_centered <- cbind(lpsa = Prostate$lpsa, covars_centered)

# Design matrix (includes intercept automatically)
X <- model.matrix(lpsa ~ ., data = Prostate_centered)
y <- Prostate_centered$lpsa
w <- rep(1, length(y))  # Equal weights

# Display dataset summary
cat("Dataset dimensions:", nrow(Prostate_centered), "observations,", ncol(X), "coefficients\n")
#> Dataset dimensions: 97 observations, 9 coefficients
cat("Response variable: lpsa (log prostate-specific antigen)\n")
#> Response variable: lpsa (log prostate-specific antigen)
cat("Predictors:", paste(colnames(X), collapse = ", "), "\n")
#> Predictors: (Intercept), lcavol, lweight, age, lbph, svi, lcp, gleason, pgg45

Prior Specification

We use comparable priors for fair comparison:

# Prior for tauBayesW (all methods)
prior_tbw <- tauBayesW::prior(
  p = ncol(X),
  type = "univariate",
  beta_mean = rep(0, ncol(X)),
  beta_cov = diag(1e6, ncol(X)),  # Very vague prior
  sigma_shape = 0.001,            # For ALD method
  sigma_rate = 0.001,
  names = colnames(X)
)

# Prior for bayesQR 
prior_bqr <- bayesQR::prior(
  lpsa ~ ., data = Prostate_centered,
  beta0 = rep(0, ncol(X)), 
  V0 = diag(1e6, ncol(X))  # Same vague prior
)

print("Prior specifications:")
#> [1] "Prior specifications:"
print(prior_tbw)
#> $b0
#> (Intercept)      lcavol     lweight         age        lbph         svi 
#>           0           0           0           0           0           0 
#>         lcp     gleason       pgg45 
#>           0           0           0 
#> 
#> $B0
#>             (Intercept) lcavol lweight   age  lbph   svi   lcp gleason pgg45
#> (Intercept)       1e+06  0e+00   0e+00 0e+00 0e+00 0e+00 0e+00   0e+00 0e+00
#> lcavol            0e+00  1e+06   0e+00 0e+00 0e+00 0e+00 0e+00   0e+00 0e+00
#> lweight           0e+00  0e+00   1e+06 0e+00 0e+00 0e+00 0e+00   0e+00 0e+00
#> age               0e+00  0e+00   0e+00 1e+06 0e+00 0e+00 0e+00   0e+00 0e+00
#> lbph              0e+00  0e+00   0e+00 0e+00 1e+06 0e+00 0e+00   0e+00 0e+00
#> svi               0e+00  0e+00   0e+00 0e+00 0e+00 1e+06 0e+00   0e+00 0e+00
#> lcp               0e+00  0e+00   0e+00 0e+00 0e+00 0e+00 1e+06   0e+00 0e+00
#> gleason           0e+00  0e+00   0e+00 0e+00 0e+00 0e+00 0e+00   1e+06 0e+00
#> pgg45             0e+00  0e+00   0e+00 0e+00 0e+00 0e+00 0e+00   0e+00 1e+06
#> 
#> $c0
#> [1] 0.001
#> 
#> $C0
#> [1] 0.001
#> 
#> attr(,"class")
#> [1] "bqr_prior"
print(prior_bqr)
#> $beta0
#> [1] 0 0 0 0 0 0 0 0 0
#> 
#> $V0
#>        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9]
#>  [1,] 1e+06 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
#>  [2,] 0e+00 1e+06 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
#>  [3,] 0e+00 0e+00 1e+06 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
#>  [4,] 0e+00 0e+00 0e+00 1e+06 0e+00 0e+00 0e+00 0e+00 0e+00
#>  [5,] 0e+00 0e+00 0e+00 0e+00 1e+06 0e+00 0e+00 0e+00 0e+00
#>  [6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+06 0e+00 0e+00 0e+00
#>  [7,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+06 0e+00 0e+00
#>  [8,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+06 0e+00
#>  [9,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+06
#> 
#> $shape0
#> [1] 0.01
#> 
#> $scale0
#> [1] 0.01
#> 
#> $method
#> [1] "QRc"
#> 
#> attr(,"class")
#> [1] "bayesQR.prior"

Model Fitting

bayesQR (Reference Method)

cat("Fitting bayesQR model...\n")
#> Fitting bayesQR model...
set.seed(12345)

fit_bqr <- bayesQR(
  lpsa ~ ., 
  data = Prostate_centered,
  quantile = 0.75, 
  ndraw = 60000, 
  prior = prior_bqr, 
  keep = 5  # Thin by 5
)
#> Current iteration :
#> [1] 500
#> Current iteration :
#> [1] 1000
#> Current iteration :
#> [1] 1500
#> Current iteration :
#> [1] 2000
#> Current iteration :
#> [1] 2500
#> Current iteration :
#> [1] 3000
#> Current iteration :
#> [1] 3500
#> Current iteration :
#> [1] 4000
#> Current iteration :
#> [1] 4500
#> Current iteration :
#> [1] 5000
#> Current iteration :
#> [1] 5500
#> Current iteration :
#> [1] 6000
#> Current iteration :
#> [1] 6500
#> Current iteration :
#> [1] 7000
#> Current iteration :
#> [1] 7500
#> Current iteration :
#> [1] 8000
#> Current iteration :
#> [1] 8500
#> Current iteration :
#> [1] 9000
#> Current iteration :
#> [1] 9500
#> Current iteration :
#> [1] 10000
#> Current iteration :
#> [1] 10500
#> Current iteration :
#> [1] 11000
#> Current iteration :
#> [1] 11500
#> Current iteration :
#> [1] 12000
#> Current iteration :
#> [1] 12500
#> Current iteration :
#> [1] 13000
#> Current iteration :
#> [1] 13500
#> Current iteration :
#> [1] 14000
#> Current iteration :
#> [1] 14500
#> Current iteration :
#> [1] 15000
#> Current iteration :
#> [1] 15500
#> Current iteration :
#> [1] 16000
#> Current iteration :
#> [1] 16500
#> Current iteration :
#> [1] 17000
#> Current iteration :
#> [1] 17500
#> Current iteration :
#> [1] 18000
#> Current iteration :
#> [1] 18500
#> Current iteration :
#> [1] 19000
#> Current iteration :
#> [1] 19500
#> Current iteration :
#> [1] 20000
#> Current iteration :
#> [1] 20500
#> Current iteration :
#> [1] 21000
#> Current iteration :
#> [1] 21500
#> Current iteration :
#> [1] 22000
#> Current iteration :
#> [1] 22500
#> Current iteration :
#> [1] 23000
#> Current iteration :
#> [1] 23500
#> Current iteration :
#> [1] 24000
#> Current iteration :
#> [1] 24500
#> Current iteration :
#> [1] 25000
#> Current iteration :
#> [1] 25500
#> Current iteration :
#> [1] 26000
#> Current iteration :
#> [1] 26500
#> Current iteration :
#> [1] 27000
#> Current iteration :
#> [1] 27500
#> Current iteration :
#> [1] 28000
#> Current iteration :
#> [1] 28500
#> Current iteration :
#> [1] 29000
#> Current iteration :
#> [1] 29500
#> Current iteration :
#> [1] 30000
#> Current iteration :
#> [1] 30500
#> Current iteration :
#> [1] 31000
#> Current iteration :
#> [1] 31500
#> Current iteration :
#> [1] 32000
#> Current iteration :
#> [1] 32500
#> Current iteration :
#> [1] 33000
#> Current iteration :
#> [1] 33500
#> Current iteration :
#> [1] 34000
#> Current iteration :
#> [1] 34500
#> Current iteration :
#> [1] 35000
#> Current iteration :
#> [1] 35500
#> Current iteration :
#> [1] 36000
#> Current iteration :
#> [1] 36500
#> Current iteration :
#> [1] 37000
#> Current iteration :
#> [1] 37500
#> Current iteration :
#> [1] 38000
#> Current iteration :
#> [1] 38500
#> Current iteration :
#> [1] 39000
#> Current iteration :
#> [1] 39500
#> Current iteration :
#> [1] 40000
#> Current iteration :
#> [1] 40500
#> Current iteration :
#> [1] 41000
#> Current iteration :
#> [1] 41500
#> Current iteration :
#> [1] 42000
#> Current iteration :
#> [1] 42500
#> Current iteration :
#> [1] 43000
#> Current iteration :
#> [1] 43500
#> Current iteration :
#> [1] 44000
#> Current iteration :
#> [1] 44500
#> Current iteration :
#> [1] 45000
#> Current iteration :
#> [1] 45500
#> Current iteration :
#> [1] 46000
#> Current iteration :
#> [1] 46500
#> Current iteration :
#> [1] 47000
#> Current iteration :
#> [1] 47500
#> Current iteration :
#> [1] 48000
#> Current iteration :
#> [1] 48500
#> Current iteration :
#> [1] 49000
#> Current iteration :
#> [1] 49500
#> Current iteration :
#> [1] 50000
#> Current iteration :
#> [1] 50500
#> Current iteration :
#> [1] 51000
#> Current iteration :
#> [1] 51500
#> Current iteration :
#> [1] 52000
#> Current iteration :
#> [1] 52500
#> Current iteration :
#> [1] 53000
#> Current iteration :
#> [1] 53500
#> Current iteration :
#> [1] 54000
#> Current iteration :
#> [1] 54500
#> Current iteration :
#> [1] 55000
#> Current iteration :
#> [1] 55500
#> Current iteration :
#> [1] 56000
#> Current iteration :
#> [1] 56500
#> Current iteration :
#> [1] 57000
#> Current iteration :
#> [1] 57500
#> Current iteration :
#> [1] 58000
#> Current iteration :
#> [1] 58500
#> Current iteration :
#> [1] 59000
#> Current iteration :
#> [1] 59500
#> Current iteration :
#> [1] 60000

# Extract coefficient estimates
draws_bqr <- fit_bqr[[1]][["betadraw"]]
colnames(draws_bqr) <- colnames(X)
coef_bqr <- colMeans(draws_bqr, na.rm = TRUE)

cat("bayesQR completed. Effective samples:", nrow(draws_bqr), "\n")
#> bayesQR completed. Effective samples: 12000

tauBayesW Method 1: ALD (Asymmetric Laplace Distribution)

cat("Fitting tauBayesW ALD method...\n")
#> Fitting tauBayesW ALD method...
set.seed(12345)

fit_ald <- bqr.svy(
  lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45,
  data     = Prostate_centered,
  quantile = 0.75,
  method   = "ald",
  prior    = prior_tbw,
  niter    = 50000,
  burnin   = 25000,
  thin     = 5,
  print_progress = 25000
)
#> Iteration 25000 of 50000

coef_ald <- fit_ald$beta
cat("ALD method completed.\n")
#> ALD method completed.

tauBayesW Method 2: Score Likelihood

cat("Fitting tauBayesW Score method...\n")
#> Fitting tauBayesW Score method...
set.seed(12345)

fit_score <- bqr.svy(
  lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45,
  data     = Prostate_centered,
  quantile = 0.75,
  method   = "score",
  prior    = prior_tbw,
  niter    = 50000,
  burnin   = 25000,
  thin     = 5,
  print_progress = 25000
)
#> Iteration 0 of 50000
#> Iteration 25000 of 50000

coef_score <- fit_score$beta
cat("Score method completed.\n")
#> Score method completed.

tauBayesW Method 3: Approximate Method

cat("Fitting tauBayesW Approximate method...\n")
#> Fitting tauBayesW Approximate method...
set.seed(12345)

fit_approximate <- bqr.svy(
  lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45,
  data     = Prostate_centered,
  quantile = 0.75,
  method   = "approximate",
  prior    = prior_tbw,
  niter    = 50000,
  burnin   = 25000,
  thin     = 5,
  print_progress = 25000
)
#> Iteration 0 of 50000
#> Iteration 25000 of 50000

coef_approximate <- fit_approximate$beta
cat("Approximate method completed.\n")
#> Approximate method completed.
# Combine coefficient estimates into one table
results_table <- data.frame(
  Method       = c("bayesQR", "tauBayesW - ALD", "tauBayesW - Score", "tauBayesW - Approximate"),
  Intercept    = c(coef_bqr["(Intercept)"], 
                   coef_ald["(Intercept)"], 
                   coef_score["(Intercept)"], 
                   coef_approximate["(Intercept)"]),
  lcavol       = c(coef_bqr["lcavol"], 
                   coef_ald["lcavol"], 
                   coef_score["lcavol"], 
                   coef_approximate["lcavol"]),
  lweight      = c(coef_bqr["lweight"], 
                   coef_ald["lweight"], 
                   coef_score["lweight"], 
                   coef_approximate["lweight"]),
  age          = c(coef_bqr["age"], 
                   coef_ald["age"], 
                   coef_score["age"], 
                   coef_approximate["age"]),
  lbph         = c(coef_bqr["lbph"], 
                   coef_ald["lbph"], 
                   coef_score["lbph"], 
                   coef_approximate["lbph"]),
  svi          = c(coef_bqr["svi"], 
                   coef_ald["svi"], 
                   coef_score["svi"], 
                   coef_approximate["svi"]),
  lcp          = c(coef_bqr["lcp"], 
                   coef_ald["lcp"], 
                   coef_score["lcp"], 
                   coef_approximate["lcp"]),
  gleason      = c(coef_bqr["gleason"], 
                   coef_ald["gleason"], 
                   coef_score["gleason"], 
                   coef_approximate["gleason"]),
  pgg45        = c(coef_bqr["pgg45"], 
                   coef_ald["pgg45"], 
                   coef_score["pgg45"], 
                   coef_approximate["pgg45"])
)

# Nicely formatted table
knitr::kable(results_table, digits = 3, caption = "Coefficient Estimates at quantile 0.75")
Coefficient Estimates at quantile 0.75
Method Intercept lcavol lweight age lbph svi lcp gleason pgg45
bayesQR 3.041 0.553 0.378 -0.020 0.112 0.897 -0.051 0.052 0.003
tauBayesW - ALD 2.896 0.566 0.143 -0.013 0.112 0.864 -0.066 -0.109 0.007
tauBayesW - Score 2169.720 151.538 251.179 -3.213 -238.876 -58.774 -42.352 54.445 22.435
tauBayesW - Approximate 2.930 0.569 0.082 -0.013 0.116 0.933 -0.057 -0.103 0.006

Model Summary and Print Methods

Let’s examine the output from tauBayesW methods using print() and summary():

Summary Methods (Detailed Diagnostics)

cat("=== ALD Method Summary ===\n")
#> === ALD Method Summary ===
summary(fit_ald)
#>             Length Class      Mode     
#> beta            9  -none-     numeric  
#> draws       50000  -none-     numeric  
#> accept_rate     1  -none-     numeric  
#> n_chains        1  -none-     numeric  
#> warmup          1  -none-     numeric  
#> thin            1  -none-     numeric  
#> runtime         1  -none-     numeric  
#> call           10  -none-     call     
#> method          1  -none-     character
#> quantile        1  -none-     numeric  
#> prior           4  bqr_prior  list     
#> terms           3  terms      call     
#> model           9  data.frame list     
#> formula         3  formula    call

cat("\n=== Score Method Summary ===\n")
#> 
#> === Score Method Summary ===
summary(fit_score)
#>             Length Class      Mode     
#> beta            9  -none-     numeric  
#> draws       45000  -none-     numeric  
#> accept_rate     1  -none-     numeric  
#> n_chains        1  -none-     numeric  
#> warmup          1  -none-     numeric  
#> thin            1  -none-     numeric  
#> runtime         1  -none-     numeric  
#> call           10  -none-     call     
#> method          1  -none-     character
#> quantile        1  -none-     numeric  
#> prior           4  bqr_prior  list     
#> terms           3  terms      call     
#> model           9  data.frame list     
#> formula         3  formula    call

cat("\n=== Approximate Method Summary ===\n")
#> 
#> === Approximate Method Summary ===
summary(fit_approximate)
#>             Length Class      Mode     
#> beta            9  -none-     numeric  
#> draws       45000  -none-     numeric  
#> accept_rate     1  -none-     numeric  
#> n_chains        1  -none-     numeric  
#> warmup          1  -none-     numeric  
#> thin            1  -none-     numeric  
#> runtime         1  -none-     numeric  
#> call           10  -none-     call     
#> method          1  -none-     character
#> quantile        1  -none-     numeric  
#> prior           4  bqr_prior  list     
#> terms           3  terms      call     
#> model           9  data.frame list     
#> formula         3  formula    call

Convergence Diagnostics

# Check acceptance rates and other diagnostics
cat("Model Diagnostics Summary:\n")
#> Model Diagnostics Summary:
cat("========================\n\n")
#> ========================

methods <- c("ALD", "Score", "Approximate")
fits <- list(fit_ald, fit_score, fit_approximate)

for (i in seq_along(methods)) {
  cat(sprintf("%s Method:\n", methods[i]))
  
  if (!is.null(fits[[i]]$accept_rate)) {
    cat(sprintf("  Acceptance rate: %.3f\n", fits[[i]]$accept_rate))
  }
  
  if (!is.null(fits[[i]]$draws)) {
    n_draws <- nrow(as.matrix(fits[[i]]$draws))
    cat(sprintf("  Posterior samples: %d\n", n_draws))
  }
  
  if (!is.null(fits[[i]]$final_scale)) {
    cat(sprintf("  Final MCMC scale: %.4f\n", fits[[i]]$final_scale))
  }
  
  cat("\n")
}
#> ALD Method:
#>   Acceptance rate: NA
#>   Posterior samples: 5000
#> 
#> Score Method:
#>   Acceptance rate: 0.814
#>   Posterior samples: 5000
#> 
#> Approximate Method:
#>   Acceptance rate: 0.141
#>   Posterior samples: 5000